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(Dated:) 

A simple and efficient method for characterization of multidimensional Gaussian states is sug- 
gested and experimentally demonstrated. Our scheme shows analogies with tomography of finite 
dimensional quantum states, with the covariance matrix playing the role of the density matrix and 
homodyne detection providing Stern-Gerlach-like projections. The major difference stems from a 
different character of relevant noises: while the statistics of Stern-Gerlach-like measurements is gov- 
erned by binomial statistics, the detection of quadrature variances correspond to statistics. For 
Gaussian and near Gaussian states the suggested method provides, compared to standard tomog- 
raphy techniques, more stable and reliable reconstructions. In addition, by putting together recon- 
struction methods for Gaussian and arbitrary states, we obtain a tool to detect the non-Gaussian 
character of optical signals. 



I. INTRODUCTION 

Gaussian states are building blocks of quantum infor- 
mation processing with continuous variables. In fact, 
Gaussian states can be generated and processed by means 
of linear operations. On the other hand, successful im- 
plementation of quantum information protocols requires 
efficient tools for the analysis and characterization of the 
quantum states involved in the experiment [l| . The stan- 
dard approach to characterize optical continuous vari- 
ables states is to perform homodyne measurements Q 
and reconstruct the measured quantum state on a sub- 
space of the infinitely dimensional space describing a 
multi-mode quantum harmonic oscillator This ap- 
proach is rather general and, in principle, can be used 
for reconstructing both Gaussian 0, 0, M, M, ^^'^ iion- 
Gaussian [ll[ states. On the other hand, such a gen- 
eral procedure becomes inherently inefficient once the set 
of possible states is restricted. For example, the knowl- 
edge about the Gaussian character of the measured state 
contributes a lot of prior information about the measured 
subject. This information can be used for reducing the 
number of relevant unknown parameters, making the re- 
construction (data inversion) much simpler, also avoiding 
problematic issues of the standard approach, such as the 
rapid growth of reconstruction errors with the size of the 
reconstruction space [H, . 

The main idea of the present contribution is to point 
out and exploit formal analogies between the description 
of Gaussian states and that of finite-dimensional quan- 
tum states. Based on this analogy we put forward a 
simple and efficient method for reconstructing Gaussian 



states from homodyne data. Since by definition the result 
of Gaussian tomography is a Gaussian state, the qual- 
ity of Gaussian fit can also be used to assess how far is 
the measured state from the family of Gaussian states, 
thus providing an operational definition of Gaussianity 
for quantum states [l4| . This will be demonstrated upon 
the application of the method to homodyne data taken 
on states produced by an optical parametric oscillator 
close to threshold [15 1. 



The paper is structured as follows. In Section |TT] we 
review few basics facts related to Gaussian states and 
homodyne detection. This will be used in Section IIIII 
for the formulation of a reconstruction method based on 
the detection of rotated field-quadrature operator. The 
theory is then generalized to multidimensional Gaussian 
states in Section HVl Experimental results and data anal- 
ysis are reported in in Sec.|Vl whereas Section IVTl closes 
the paper with some concluding remarks. 



II. HOMODYNING GAUSSIAN STATES 



For the sake of simplicity we start by formulating the 
problem of reconstruction for single mode Gaussian field 
g, with Wigner function given by 



W{X) = ^ ' ^ '-^ , (1) 
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with X"'" = {x,y) and covariancc matrix G given by 

If 2(AX)2 {AX, AY}\ > o = i /"O 
'^-2\{AX,AY} 2{AP)^ J - 2\i 

■ (2) 

where X — {a + a^)/\/2 and Y = i{a^ — a)/\/2, arc the 
quadrature operators, a and being mode operators. If 
the homodyne detector is set to measure the quadrature 
X{9) = XcosO + y sinfl, the positive operator valued 
measure (POVM) n,,(x, 6*) associated with its reahstie 
measurement is the Gaussian convolution of the quadra- 
ture projectors, 
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where 6"^ = {1 — r])/2ri, and ry is the quantum efficiency 
of the involved (linear) photodetectors. In the Fock basis 
we have 



\y)i> = ( z 
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and the distribution of homodyne outcomes is given by 
p^ix,9)^Tr[gUr,{x,e)] (5) 
^e^p{-{x-[RoX]y/{2af,)} 
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where = r]{M22T:>ct[G]+S^), M22 = [M]22, with M 
Rg G~^R 0, and the rotation matrix is defined as 
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III. ESTIMATION OF COVARIANCE MATRIX 

In this Section, we address the problem how to esti- 
mate efficiently the covariance matrix G. For the sake of 
simplicity let us set the coherent part of the signal to zero 
and 77 = 1. Introducing a unit vector {u\ = (cos 0, sin 0) 
parametrized by the phase of the local oscillator and 
noticing that (M|G|tt) = A/22Det[G], the sampled prob- 
ability density ([6]) may be conveniently expressed in the 
form 



p{x,9) 



v/Stt {u\G\u] 



: exp 



2(u\G\u) 
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There are several striking similarities between the recon- 
struction of spin 1/2 states and covariance matrices. In- 
deed, in both cases the state is described by 2 x 2 matri- 
ces. The density matrix of the spin 1/2 state is a hermi- 
tian, unit trace, semi-positive matrix, thus leaving three 
free parameters for its full description. The covariance 



matrix is real, symmetric matrix, again fully described by 
three parameters, constrained to the relation Pro- 
jections of a spin density matrix are conveniently sam- 
pled in a Stern-Gerlach experiment. Similarly, detections 
of quadrature variables in a homodyne experiment pro- 
vide sampling of quadrature variances representing pro- 
jections of a covariance matrix. Indeed, estimated matrix 
G appears only in the variance of distribution ([7]). Single 
detected event do not say too much, but repeated detec- 
tions do. Provided that homodyne detection is repeated 
n times for the same setting, Xi being the detected re- 
sults, the likelihood of G reads 



£(G|{xJ) 
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Consider now the statistics of the random variable 



(8) 



y 



which is given by 



Paiy) = / dxi...dx„5{y 
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It is easy to see that for Gaussian states the fluctuations 
of y variable are governed by the well-known distri- 
bution 
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Upon maximizing the likelihood (|10p . the quadrature 
variance, which plays the role of a projection of the co- 
variance matrix, cr^(G) = (it|G|M), is estimated as 



a'{G) 



(11) 



This establishes a formal analogy between estimations 
of spin states and Gaussian covariance matrices. In the 
former case the probability of finding the spin "up" is 
sampled by the number of particles defiected upwards in 
a Stern-Gerlach apparatus. Similarly, for Gaussian states 
the variance of a quadrature distribution (i.e. a projec- 
tion of the covariance matrix) is sampled by properly 
normalized sum of squares of the detected quadrature 
values. 

Having mentioned similarities between the two prob- 
lems let us also identify three important differences be- 
tween the estimation of a spin 1/2 system and estimation 
of a Gaussian state: 

(i) the two problems have different underlying statics: 
binomial statistics for yes-no spin data and dis- 
tribution for sampled variances; 

(ii) there are slightly different constraints on the es- 
timated quantities: semi-positivity of density ma- 
trix for spins and uncertainty relations constraint, 
stated in the form of the matrix inequality for 
Gaussian states; 
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(iii) the measured quantities have different nature: for 
spins the sampled probabilities are expectation val- 
ues of Stern-Gerlach outcomes whereas for Gaus- 
sian states we need sampled variances of the homo- 
dyne data, i.e the information is contained in the 
data noise. 

Notice that properties (ii) and (iii) will be important for 
the discussion of our main result below. 

In order to get a unique reconstruction of the mea- 
sured Gaussian state, quadrature measurement should 
be repeated with different settings of phases 9h of 
the strong local oscillator. Let us assume that each 
quadrature X{9h) is sampled Uh times with the results 
Xh,i, Xh,2, ■ ■ ■ Xh.uh, and denote yh = Y.k '^l.k- The cor- 
responding log-likelihood reads 

log/:(G) = -i^n,loga^(G)-^— 1^, (12) 

h h h\ ) 

where ^^(G) = {uh\G\uh). A reconstruction of the sig- 
nal covariance matrix is then obtained by maximizing 
the likelihood hmction subject to the constraint G > fZ. 
Let us first consider the extremal equation for the covari- 
ance matrix. It is easily found by taking a derivative of 
log£(G) with respect to matrix G. This yields 

RG = DC, (13) 

where 

^-i:(i^:^K)Ki, (14) 

As there is again a direct analogy between the extremal 
equation for G and the corresponding extremal equation 
for the maximum-likely spin state, methods developed for 
generic quantum estimation can be employed, with some 
caution, to solve the extremal equation (fT3| by iterations 
p^ . A form suitable to iterations can be foimd noting 
that by Hermicity of G, R and D both G = D^^RG and 
G = GRD^^ simultaneously hold and may be combined 
into a single extremal matrix equation 

G = D-^RGRD~\ (16) 

which is our main formal result. This matrix equation 
may be iterated starting e.g. with the covariance matrix 
of the vacuum state. In addition the form of Eq. p6|) 
guarantees that semi-positivity, G > 0, of the covariance 
matrix is preserved after each iteration step. 

Some discussion of this algorithm is now in order. As 
we have already mentioned above condition G > would 
be typical for density matrices. Covariance matrices obey 
a more complicated relation G > fl. Though in princi- 
ple the algorithm (|16p may converge to a covariance ma- 
trix violating Hcisenberg uncertainty relations, we will 



present physical arguments showing that in practice this 
almost never happens. 

First, let us discuss the estimation of a spin state from 
noisy data. Ignoring for a while the constraint p > 0, the 
probability that the best fit of measured data is provided 
by a non physical matrix depends on the distance of p 
from the boundary of the convex set of density matri- 
ces and also on the character of the noisy data. When 
the true state lies close to the boundary and/or the mea- 
sured relative frequencies differ a lot from the theoretical 
probabilities the chances are high that the data are best 
explained by a "density matrix" having at least one nega- 
tive eigenvalue. In addition, states close to the boundary 
are easy to produce e.g. by a projection. Hence the 
condition p > is an essential part of the spin state to- 
mography, which must be incorporated in any meaningful 
state reconstruction protocol. 

The situation is quite different in Gaussian state to- 
mography where the boundary of physical covariance 
matrices is made by minimum uncertainty states, which 
rarely appear in practice. In fact, in order to generate a 
minimum uncertainty state one has to eliminate all the 
sources of noise from the preparation procedure. In addi- 
tion, variances (noise) rather then probabilities are fitted 
from data in Gaussian tomography. Hence, any addi- 
tional noise, being of statistical origin or due to some 
imperfection in the measurement, will have the effect of 
moving the reconstructed state further away from the 
boundary. For these reasons the condition G > H is 
not likely to be violated by G reconstructed from real 
data. This justifies the use of the simple reconstruction 
algorithm of Eq. pG)) . Cast in another way, a possible 
violation of the condition G > H can be seen as a chal- 
lenge for experimenters. Indeed, to obtain such a data 
set, one must be able to prepare and control quantum 
states lying very close to the minimum uncertainty states, 
and to measure them with sufficient precision. This re- 
quirement is obviously much more demanding then, for 
example, the generation and observation of squeezing in 
one quadrature without paying attention to the excess 
noise introduced in the conjugated observable. In case 
of such an advanced experiment the condition G > 
can taken into account by means of a simple two-step 
strategy. Provided G > fl the reconstructed covariance 
matrix is physically sound and the search is over. In the 
opposite case a search for the maximum-likely covariance 
matrix should be repeated among the set of minimal un- 
certainty states. 



IV. RECONSTRUCTION OF MULTI-MODE 
GAUSSIAN STATES 

The reconstruction method illustrated in the previous 
Section may be straightforwardly generalized to the prob- 
lem of estimating the covariance matrix of multi-mode 
Gaussian states. The strategy is the following. The 
modes will be mixed by means of controlled beam split- 
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ter and the homodyne detection will be done on the se- 
lected output mode. Provided that covariance matrix has 
been sampled by sufficient number of projections, it can 
be reconstructed. In order to apply the ML procedure 
developed above it is enough to show how variances of 
the detected modes can be derived from the generic co- 
variance matrix. Without loss of generality the method 
will be illustrated on the example of two-mode case since 
the extension of the model to higher dimensional case is 
st r aight forward . 

Let Gin be the 4x4 covariance matrix of a two-mode 
Gaussian state with entries 



1 



Gki = = -{{Qk,Qi}) - {Qi){Qk) 



(17) 



where Q = (Xi, X2, ^2), {Qk,Qi} = QkQi + QiQk, 
and = -i=(afc + aj.) and Yfe = '[^("'k O'l) are the 
quadrature operators of mode k. Moreover, one should 
have 



Gin >ft®n 



n 








ft 



(18) 



being the 2x2 null matrix. If the two modes are 
mixed at a beam splitter with transmissivity r = cos^ 1?, 
then the emerging two-mode state is still Gaussian and 
its covariance matrix is given by: 



Gi, 



Go 



where 



SbsW 



cos l9 I2 


sin-d I2 


-sint?l2 


cos-i? I2 



(19) 



(20) 



is the symplcctic transformation associated with the 
beam splitter and I2 is the 2x2 identity matrix. In the 
case of a balanced beam splitter (i? = 7r/4) the two out- 
put modes, say 61 and 62, correspond bi = (ai + 02) /V^ 
and 62 = (oi — 02) /V^, respectively. In this case we can 
write the evolved covariance matrix as follows: 



Gn 





C 







with 
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G24 



[Gi + G2 ± (G3 + G4)] , 
[G2 - Gi + (G3 - G4)] , 



G2 = 



G33 G34 
G34 G44 



G4 = G 
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(23) 

(24) 
(25) 



and the elements Gki are defined by Eq. (fT7|) . Ma- 
trices Gi and G2 correspond to the covariance matri- 
ces of the reduced single- mode input states, i.e., when 
the other one is neglected. We have chosen the no- 
tation in such a way that ^4^+^ and A^~'> are the co- 
variance matrices of the single-mode states correspond- 
ing to modes bi and 62, respectively. Homodyne de- 
tection on mode bi can be simply described by ex- 
tension of vector \u) into higher dimension as {{u\ = 
(cos0, sin0, 0, 0); similarly, in the case of mode 62 one 
has to use the vector {{v\ = (0, 0, cosS, sin^^). Hence the 
detected variance can be easily represented by expres- 
sion {{w\Sbs{'>^) GinS^g{-d) \w)) , with w = u,v. Con- 
sequently, the reconstruction can be done by means of 
the same reconstruction algorithm as above considering 
index h as multi-index for phase of homodyne detection 
6k and the beam splitter transmissivity t = cos^ z?: 



SEsWluk)). 



(26) 



In particular, two configuration of beam splitter are suf- 
ficient for successful reconstruction of 2-mode covariance 
matrix, namely the detection on the free (uncoupled) 
modes with S'bs(O) = 1 [corresponding to covariance 
matrices ([24]) ] and mixing with symmetric beam split- 
ter [corresponding to covariance matrices (HH)]. Taking 
into account the efficiency, the extremal equation for co- 
variance matrix reads 



D{w)G = R{w)G, w = u,v 



(27) 



where 



and 



Riw)=J2 



{{wk\SBsWGS^smwk))+6^^ 

X S]^smwk)){{wk\SBsW (28) 



-^[{{wk\SBsWGs^^m'^k))+sf^y 

xS^sW\'^k)){{wk\SBs{^). 
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V. EXPERIMENTS AND DATA ANALYSIS 

The method exposed in the previous section has been 
applied to homodyne data samples acquired for two dif- 
ferent kind of input state: a (definitely Gaussian) refer- 
ence vacuum state and slightly non Gaussian state, close 
to squeezed thermal vacuum state, generated by an I- 
type optical parametric oscillator (OPO) that operates 
close to the threshold. In the following, we give a brief 
description of the experimental setup [l5| with particular 
attention payed to the homodyne detector and the data 
acquisition system. The setup can be divided into three 
distinct blocks: the state source, a below threshold fully 
degenerated OPO, the detector, a quantum homodyne 
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Homodyne 




FIG. 1: Experimental setup. 



and related electronics, and the, PC based, acquisition 
board (see Fig. [1]). 

The source is a below threshold fully degenerate OPO 
based on a type-I LlNbOsiMgO crystal and pumped at 
532 nm [isj . The quantum homodyne detector, whose 
reliability has been proved in different experiments 0, 
[T5| . shows an overall quantum efficiency of r/ = 0.88 ± 
0.02. To avoid any influence on the statistics of the data, 
the electronic noise floor is kept w 15 dB below the shot 
noise. Mode matching at the BS has been accomplished 
by spatially flltering the local oscillator (LO) beam by a 
mode cleaner cavity whose geometrical properties match 
the OPO's ones. The relative phase between the LO 
and the beam exiting the OPO is scanned by a piezo 
mounted mirror to which a linear voltage ramp is applied; 
the ramp is adjusted so to have a 2it phase span in the 
acquisition time. The 27r variation necessary for a full 
state reconstruction is selected in the central region of 
the piezo range so to optimize the linear response of the 
piezo stack. 

Data sampling is moved away from the laser carrier 
frequency by mixing the homodyne current with a sinu- 
soidal signal of frequency = 3 MHz. The resulting cur- 
rent is low-pass filtered with a bandwidth B = 1 MHz, 
and eventually sampled by a digital acquisition PC based 
module (Gage 14100) able to acquire up to IM-points 
per run with 14 bits resolution. For each state, 1 048 308 
quadratures values were measured. 

Measured quadrature values normalized with respect 
to vacuum fluctuations are shown in Fig. [2] To facili- 
tate our analysis the data were grouped into 31 phase 
bins prior to reconstruction. The purpose of the experi- 
ment was to perform Gaussian tomography and compare 
the results with that of the standard homodyne tomog- 
raphy p^j where no assumption about the Gaussianity 
of the measured state were made. The comparison of 
the reconstructed states may be seen as a test of Gaus- 
sian behavior. Such an approach is superior to a test 
based on evaluating certain quadrature moments due to 
the fact that all the relevant moments are encoded in the 
reconstructed state. 

The Gaussian analysis starts with evaluating quanti- 
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FIG. 2: Raw homodyne data for the vacuum state (top) and 
for the OPO output close to threshold (bottom). 

ties Dh = X]fe ^ft, j foi' s^ch phase bin ft, = 1, 31. These 
are used to construct operator R appearing in the ex- 
tremal equation for covariance matrix Eq. (|13p . The 
following covariance matrices were found for measured 
vacuum (V) and OPO (O) data: 

/ 0.50 0.00 \ / 2.38 -0.53 \ 

y 0.00 0.50 / \ -0.53 0.55 j ^ ' 

Matrix Gq indicates that the reconstructed Wigner func- 
tion of the OPO state is slightly rotated with respect to 
the x-p axes, which are defined by the the phase of the 
local oscillator. The squeezed nature of the input signal 
is revealed by diagonalizing matrix Go- With that vari- 
ances of 0.40 and 2.53 are found for the squeezed and 
anti-squeezed quadratures, respectively. Notice that the 
reconstructed state is not a minimum-uncertainty state, 
DetGo = 1.02 ^1/4 and so the reconstructed covari- 
ance matrix is located well inside the set of physical co- 
variance matrices, Go > To summarize the Gaus- 
sian analysis, a squeezed Gaussian state is inferred from 
experimental data whose anti-squeezed quadrature has 
been contaminated by some classical excess noise. Now 
it is interesting to compare the result of Gaussian analy- 
sis with that of the generic quantum tomography applied 
to the same data set. In this case, no assumption is made 
about the Gaussian character of the measured signal and 
the reconstruction is done on an iV-dimensional trun- 
cated Fock space using the maximum-likelihood tech- 
nique [TtI . Data have been grouped into 31 phase bins. 
In addition, detections within each phase bin have been 
divided into 31 quadrature bins. The set of quadrature 
projectors corresponding to the 31 x 31 = 961 pairs of 
phase and quadrature values comprises a POVM describ- 
ing the tomography scheme. Prior to reconstruction, the 
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dimension N of the reconstruction space must be cho- 
sen. Notice that this parameter has no analogy in the 
Gaussian analysis above. 




FIG. 3: Hilbert-Schmidt distances of state reconstructions 
done in Ai'-dimensional Fock spaces with respect to the refer- 
ence state reconstructed in dimension N = 30. 

In order to see the effects of the truncation on the re- 
constructed state, the procedure has been repeated for 
different values of A^. The convergence of the generic re- 
construction with growing reconstruction space is shown 
in Fig. [3l where the Hilbert-Schmidt distance dns = 
Tr{(p7v — P3o)^} between a reconstruction from a trun- 
cated space and a reference reconstruction from a large 
space is plotted function of N. This figure indi- 
cates that from = 10 on the reconstructed state does 
not change significantly when the reconstruction space is 
enlarged. This conclusion, however, turns out to be in- 
correct when one is interested in qualitative, rather than 
quantitative properties of the reconstructed state. For 
instance, the shape of the reconstructed Wigner func- 
tion may help to classify the signal as being Gaussian or 
non-Gaussian. Similarly, the presence of negative values 
of W may be taken as a sign of non-classical behavior. 
In turn, making qualitative statements from quantitative 
results of tomography may be delicate and challenging. 
This is illustrated in Fig. |4l where we report the Wigner 
function and the corresponding y-cut for different val- 
ues of the truncation dimension for the reconstruction in 
the Fock space. Notice that what may seem as a sign 
of non-classical interference in Fig. SJ^a) is actually an 
artifact due to small value of the reconstruction dimen- 
sion. Remarkably, this happens even though the inspec- 
tion of the elements of pi5 and p25 in the Fock basis 
shows negligible differences (see also Fig. [S]), since these 
tiny differences get amplified when switching to Wigner 
representation. Also plots in Fig. [DJb) distinctly shows 
non-Gaussian shapes and the artifacts are present un- 
til dimension of TV w 25 is reached. In this region the 
generic reconstruction starts to be consistent with the 
Gaussian results yielding very similar values of quadra- 
ture variances. More than 600 free parameters must be 
introduced in the standard tomography to achieve the 
quality of the simple Gaussian fit of OPO data! 

There remains a question whether the deviations from 



Gaussian character that may be observed in the generic 
reconstruction even for large dimensions N are real fea- 
tures of the measured state or not. This possibility may 
be confirmed or discarded by analyzing the scores (likeli- 
hoods) achieved by the Gaussian and generic reconstruc- 
tions respectively. Here, a direct comparison is not possi- 
ble due to different natures of involved data. Instead, one 
may determine how much the quality of the fit degrades 
by replacing vacuum data with the more noisy OPO sig- 
nal. For Gaussian analysis we find that the log-likelihood 
of the OPO reconstruction is smaller by a factor 3.66 
compared to the vacuum reconstruction. Standard to- 
mography yields a smaller factor of 3.16, which indicates 
that a better fit of the OPO signal generated near the 
threshold is obtained by abandoning the Gaussian hy- 
pothesis. 

In order to test the null hypothesis for our samples 
(i.e., data normally distributed) we grouped homodyne 
data into phase bins made of 10 000 outcomes each and 
applied "normality tests" to both the samples of Fig s. [H 
In particular, we have used the Jarque-Bera (JB) [l9[ and 
the Shapiro- Wilk (SW) normality tests The JB test 
is based on both the skewness S and the kurtosis (excess) 
K of the sample, i.e 



S 



(M2)3/2 



K 



3, 



(31) 



where fi^ = jf 'I2h=ii-'^h — x)'' is the fc-tli central moment, 
N is the number of data and x their mean. The JB 
statistics is then defined as: 
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(32) 



Since Wjb has an asymptotic distribution with two 
degrees of freedom [l^l, one can use it to test the null 
hypothesis of normality. In particular, by choosing a 
significance level 0.05, one rejects the null hypothesis if 
Wm > 5.99. The test statistics Wsw of the SW test 
reads as follows: 



=1 ahX(h) 



(33) 



where ccj/j-j are the ordered sample values (x(/j) is the /i-th 
smallest value) and ah are constants generated from the 
means, variances and covariances of the order statistics of 
a sample of size N from a normal distribution . In this 
case one rejects the null hypothesis of normality within 
a significance interval of 0.05, if p-Wsw < 0.05, where 
the p-Wsw is the p- value of Wsw^ i-e- the probability of 
obtaining a result at least as extreme as the one that was 
actually observed, given that the Gaussian hypothesis is 
true. The two test statistics are shown in Fig. [51 where 
we plotted the variance, Wjb andp-VFsw associated with 
the data of Figs. (5] As one can see, the vacuum data 
arc normally distributed, whereas those coming from the 
OPO close to threshold show clear deviations from the 
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FIG. 4: Wigner functions (upper row) and cuts along the y axis (bottom row) of the OPO state reconstructed in Fock spaces 
of dimensions A'' = 8 (a), = 15 (b), and A'' = 25 (c) respectively. 



Gaussian behavior (see also [3|), thus justifying the re- 
jection of Gaussian hypothesis. 



VI. CONCLUSIONS 

In this paper, upon exploiting the formal analogies 
between the description of Gaussian states and that of 
finite-dimensional quantum states, we have proposed and 
demonstrated a simple and efficient method for the re- 
construction of Gaussian states. The method have been 
tested numerically and applied to the reconstruction of 
the quantum state of the signal generated by an opti- 
cal parametric oscillator. As a matter of fact, under the 
Gaussian hypothesis there is a drastic reduction in the 
number of parameters needed to characterize a quantum 
state, and this leads to a robust reconstruction technique. 
Indeed, our results indicate that on Gaussian and near 
Gaussian states the proposed method provide a more sta- 
ble and reliable reconstructions than standard tomogra- 
phy techniques. We have also shown that putting to- 
gether results coming from the Gaussian method with 
those coming from standard tomography, we obtain a re- 
liable tool to detect the non-Gaussian character of optical 
signals. We conclude that the present method provides a 
reliable and robust tool for the characterization of quan- 
tum states which may be of interest for a rather broad 
community working in continuous-variable quantum in- 
formation processing and quantum state reconstruction. 
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FIG. 5: (Color online) Variance (circles), Wm (triangles) and 
the p-value of Wsw (p-Wsw, squares) associated with the 
homodyne data of Figs. [21 Top: vacuum data. Bottom: OPO 
close to threshold data (here p-Wsw has been magnified by a 
factor xl.5). The dashed line is the threshold value 0.05 for 
p-Wsw (xl.5 in the bottom plot); the dot-dashed hue is the 
threshold value 5.99 for Wjb- If Wjb > 5.99 and/or p-Wsw < 
0.05, then the null hypothesis (Gaussian states) is rejected. 
See the text for details. Results indicate that vacuum data 
(top) are consistent with the Gaussian hypothesis whereas for 
data coming from the OPO close to threshold the Gaussian 
hypothesis should be rejected. 
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